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We instigate the properties of the threshold contact process (TCP), a process showing an 
absorbing-state phase transition with infinitely many absorbing states, on random complex net- 
works. The finite size scaling exponents characterizing the transition are obtained in a heterogeneous 
mean field (HMF) approximation and compared with extensive simulations, particularly in the case 
of heterogeneous scale-free networks. We observe that the TCP exhibits the same critical properties 
as the contact process (CP), which undergoes an absorbing-state phase transition to a single absorb- 
ing state. The accordance among the critical exponents of different models and networks leads to 
conjecture that the critical behavior of the contact process in a HMF theory is a universal feature of 
absorbing state phase transitions in complex networks, depending only on the locality of the inter- 
actions and independent of the number of absorbing states. The conditions for the applicability of 
the conjecture are discussed considering a parallel with the susceptible-infected-susceptible epidemic 
spreading model, which in fact belongs to a different universality class in complex networks. 
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I. INTRODUCTION 

The role of the effects of a disordered substrate on 
the behavior of dynamical processes has attracted the 
attention of the statistical physics community since long 
[TJ [2J . The interest in this issue has been enhanced by 
the recent realization that the substrate of many rele- 
vant dynamical processes can be represented in terms of 
a complex network j3H5J. This realization has led to a 
renewed scientific effort emphasizing the effects of the 
disordered topology inherent to the underlying network 
[3 [7]. In this context, it has been shown that highly het- 
erogeneous contact patterns, characterized by a degree 
distribution P(k) (probability that an element or vertex 
in the network is connected to other k vertices) show- 
ing a heavy-tailed form can have a very strong impact in 
processes such as epidemic spreading [8HTD]. percolation 
[TTT [T2] , or synchronization [TSJ . 

Among the different dynamical processes considered 
in networks, an important role has been played by the 
class of processes with absorbing states [HI H5J. Ab- 
sorbing states are configurations of the system in which 
the dynamics becomes trapped forever and that usually 
imply the presence of absorbing state phase transitions 
(APTs) between an active phase and the absorbing con- 
figurations. The simplest model exhibiting an APT is 
the contact process (CP) [T5], where particles lying on a 
lattice or network undergo spontaneous annihilation and 
catalytic creation. A creation event involves a pair of 
empty and occupied nodes and is controlled by a rate A. 
The CP undergoes a continuous APT at a critical point 
A c , separating an absorbing phase with all sites empty 
from an active one with a nonzero order parameter </>, 
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defined as the average density of particles at the steady 
state. The ensuing APT is characterized in terms of a 
set of critical exponents, that define the so-called directed 
percolation (DP) universality class [T5J. The theoretical 
analysis of the CP in networks [TTl - fPT)] , based in the het- 
erogeneous mean-field (HMF) theory 0(7], has shown the 
presence of an APT in which the scaling of the relevant 
critical quantities with network size presents very strong 
corrections to scaling in scale-free (SF) networks with an 
heterogeneous degree distribution given by a power law 
form, P(k) ~ fc -7 [3J. These corrections to scaling, albeit 
very large, have been observed in large-scale numerical 
simulations [2D] . 

While simple models such as the CP are character- 
ized by a unique absorbing state devoid of all particles, 
other relevant systems present multiple absorbing con- 
figurations [21 . Examples of those range from models of 
forest fires [22J to self-organized critical sandpiles [23] . A 
basic model presenting an APT with multiple absorbing 
states is the pair contact process (PCP), where a pair 
of occupied neighbors annihilates with probability p and 
creates a new particle in one nearest neighbor (NN) of 
the pair with probability 1 — p |24j . Any configuration 
with only isolated particles is an absorbing state, which 
implies that the PCP has infinitely many absorbing con- 
figurations in the thermodynamic limit. On lattice, the 
PCP belongs to the DP universality class [15] . 

In the field of complex networks, dynamical systems 
with many absorbing states have been used to inves- 
tigate self-organized criticality and avalanches [251127] . 
while the analysis of the ensuing APT is limited to a 
few works [28] [29] . The HMF analysis of the forest fire 
model in heterogeneous networks yielded critical expo- 
nents equal to those obtained for the susceptible-infected- 
susceptible epidemic model in the same approach [28j . A 
finite-size scaling [30] analysis of the PCP was performed 
for homogeneous networks [29] . but using very limited 
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sizes, which do not allow to extract reliable conclusions 
about the model's universality class. 

In this work, we present and analyze a model pos- 
sessing an APT with many absorbing states, which is 
amenable to analytical calculations based on the HMF 
formalism and large-scale numerical simulations. The 
model, that we call the threshold contact process (TCP), 
is inspired in the PCP, with the peculiarity that the 
fermionic constraint of single occupation of vertices is 
relaxed to allow double occupation. This modification 
hugely reduces the computer algorithmic complexity and 
allows us to perform a standard HMF analysis [51 [71 131] . 
We show that the behavior of the TCP at the HMF level 
is the same of the CP, a fact that we confirm by means 
of extensive numerical simulations based in the quasi- 
stationary state method [THl HOI 122] • We thus conclude 
that, despite the non-fermionic nature of the TCP, it ex- 
tends to complex networks the fundamental inclusion of 
the PCP in the universality class of the contact process. 

We have organized our paper as follows: In Sec. [H] we 
present a brief overview of the properties of the CP on 
complex networks. Section [lIT| is devoted to describe the 
TCP model. A HMF theory and a homogeneous pair- 



approximation for TCP are developed in Sec. IV The 



simulation methods and numerical procedures used are 
described in Sec. |Vj where we also present a check of 
the HMF theory by means of numerical simulations on 
annealed networks. Simulations on real quenched net- 
works are presented in Sec. |VI| where we compare them 
with the HMF theory predictions. Finally, we draw our 
concluding remarks in Sec. |VII| 



II. THE CONTACT PROCESS ON NETWORKS 

The CP dynamics on networks is defined as follows 33 1 : 
Every node can be either occupied or empty. Occupied 
nodes become empty at a rate 1, while empty sites be- 
come occupied at rate Xn/k, where k is the number of 
NNs of the node (the node's degree), and n the number 
of occupied nodes among the NNs. 

The properties of the CP dynamics on heterogeneous 
networks have been analytically investigated within the 
framework of the heterogeneous mean-field (HMF) the- 
ory [TTHT§1 [33] , an extension of the standard mean-field 
theory taking into account the network heterogeneity. 
The HMF theory assumes that the vertex degree is the 
only relevant quantity, substitutes the topological struc- 
ture of the network by an annealed form [H] [7] , and ne- 
glects all dynamical correlations between vertices. In 
HMF theory, focus is placed on the partial densities of 
particles in vertices of degree k, denoted by 4>k, which 
satisfy the rate equation [33] 



dt 



\k[i - 4> k ] 



k> 



(1) 



a vertex of degree k is connected to a vertex of degree 
k' [34]. The solution of this equation in the station- 
ary limit leads to an APT located at the critical point 
A c = 1, irrespective of the degree correlations in the 
network. For degree uncorrelated networks, in which 
P(k'\k) = k! 'P(k')/ '(k) [4], a complete solution can be 
worked out [33j . giving an order parameter, defined as 
the average particle density, that close to the critical 
point scales as <j) = ~ (A — X C Y , a char- 

acteristic time scale r ~ (A — X c )~ v and a decay of the 
particle density at the critical point </> c (i) ~ t~ s . In SF 
networks with degree distribution P(k) ~ k^ 1 , one has 
f3 = 5 = max{l/(7 — 2), 1}, while v = 1 for any 7. 

These critical exponents, depending on the degree dis- 
tribution, are valid in the thermodynamic limit of an in- 
finite network. In finite systems, however, the transition 
to the absorbing state is better characterized in terms of a 
finite-size scaling (FSS) analysis [30]. In Refs. [T7 | H"9 l 135] 
such theory has been developed, exploiting a mapping of 
the CP in the low density regime to a one-step process. 
Within this approximation, one can build the equation 
for the probability distribution P n of observing n active 
particles in the steady state, which in the thermodynamic 
limits takes a scaling form given by [19] 



P, 



f 



(2) 



where / is a scaling function normalized as f(x)dx = 
1, fi = N/g and g = (k 2 )/(k) 2 . The critical density 
and characteristic time follow directly from this distribu- 
tion |19| and are given by 



4> ~ (Ng) 



-1/2 



and 



(N/g) 



1/2 



(3) 



For a power law degree-distribution, the factor g diverges 
with the degree cutoff k c (the largest degree present in 
the network) as g ~ k% 1 for 7 < 3 and becomes constant 
(logarithmic) for 7 > 3 (7 = 3) [36]. For a general cutoff 
scaling as k c ~ N 1 ^, the scaling laws 



0~ N~ 



and 



N a 



(4) 



ensue in the limit N — > 00, with associated critical expo- 
nents [19] 

v — -+max ( - - , ) , a = -—max | - - , ) . (5) 



2uj 



2uj 



This exponents, which define what we can call the CP 
universality class at the HMF level (the CP- HMF class), 
have been proved to provide a correct description of the 
CP in networks by means large-scale numerical simula- 
tions [20] , corroborating thus the validity of the descrip- 
tion of the CP on networks in terms of the HMF theory. 



III. THE THRESHOLD CONTACT PROCESS 



The quantity P(k'\k) is a measure of the degree corre- 
lations [1], defined as the conditional probability that 



One of the simplest models exhibiting a transition to 
infinitely many absorbing states is the PCP [24] , An 
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optimized computer simulation of PCP near the critical 
point (low density of active sites) requires the mainte- 
nance of an updated list containing the coordinates of 
all active pairs. This task becomes computationally in- 
efficient when the model is implemented on an arbitrary 
graph, more specifically during the annihilation events 
when the involved pair and all the corresponding active 
links must be removed from the list. Kockelkoren and 
Chate [37] had already realized that the fermionic con- 
straint can be counterproductive in lattices due to algo- 
rithmic complexity and, more importantly, due to strong 
deviations to scaling originated from this restriction. For 
this reason, we consider a modification of the PCP, that 
we call the threshold contact process (TCP), that is de- 
fined by the following rules: Nodes can be occupied by 0, 
1 or 2 particles (i.e. they can have an occupation number 
Cj = 0, 1, 2). Particles in a doubly occupied node can 
be both annihilated at rate 1 or either create an offspring 
particle in a randomly chosen NN at rate A. The creation 
event takes actually place only if the selected NN is empty 
or singly occupied. Because empty and singly occupied 
nodes do not react, any configuration devoid from dou- 
bly occupied sites is an absorbing state. Singly occupied 
sites do however play a fundamental role, since they form 
a backbone of connected vertices through which activity 
can spread. The TCP can be easily generalized to higher 
thresholds (maximum occupation number), h > 2. 

The computer implementation of TCP in an arbitrary 
graph can be done with the following scheme: At each 
time step, an active (doubly occupied) vertex j is ran- 
domly chosen and time is updated ast->i + At, where 
At = 1/[(1 + X)n(t)] and n(t) is the number of active 
vertices at time t. With probability p = 1/(1 + A), both 
particles of the selected vertex are annihilated, producing 
an empty vertex. With the complementary probability 
1— p = A/(l + A), one of the k neighbors of j is randomly 
selected and, if empty or singly occupied, receives a new 
particle; otherwise nothing happens and the simulation 
proceeds to the next time step. 



IV. MEAN-FIELD ANALYSIS 
A. Homogeneous mean-field theory 

In homogeneous substrates such as a regular lattice, 
all nodes are equivalent and thus the equations of mo- 
tion are the same for all of them. Letting P(oi, • • • ,<ti) 
be the probability that a cluster of I vertices has the 
configuration {<j\ ■ ■ ■ a{\, we can find the following exact 
equations of motion for single sites in the TCP dynamics: 



dP(2) 

dt 
dP(l) 

dt 



= -P(2) + AP(1,2) 
= A[P(0,2)-P(1,2)]. 



(6) 
(7) 



The equations for single nodes depend on pairs, those 
for pairs depend on triplets (as we will see below), and 



so forth. An approximate solution can be obtained by 
breaking this hierarchy of equations by means of some 
closure condition. Thus, in the simplest one-site approx- 
imation, the correlation between pairs is neglected, re- 
sulting in P(a,a') « P{a)P{a') U\. Equations Q and 
([7]) in this one-site approximation become 



dt 
dp 

~dt 



-(/) + \pcf>, 
\<t>{l-ct>-2p), 



(8) 
(9) 



where we have defined <f> = P(2), p = P(l), and obvi- 
ously P(0) = l — p — 0. In the stationary regime, Eqs. ^ 
and ([9]) yield p = 1/A and <fi = (A - 2)/A. We can there- 
fore identify the critical point A c = 2 and the DP mean 
field exponent j3 = 1 [14] . One can easily verify that the 
characteristic time r scales with DP mean field exponent, 
r ~ (A — Ac) -1 , while the density of active sites decays 
with time at criticality as 4> c ~ t . 

For the case of the CP, it has been shown that a ho- 
mogeneous pair approximation (HPA) [38] . considering 
the hierarchy of equations up to second order, yields bet- 
ter estimates of the critical point even in the realistic 
case of heterogeneous networks [2UJ We can also 

follow this HPA approach for the TCP. Firstly, we ob- 
serve that in the stationary regime, Eqs. ([6| and ^ give 
P(2, 0) = P(2, 1) = 4>/X. Moreover, using the normal- 
ization condition P(0, 2) + P(l, 2) + P(2, 2) = P(2), we 
have P(2, 2) = 1 — 2(j)/X. The equations for the remaining 
independent pairs are 

dP(o, o) = 2p ^ Q) _ 2 q_^2 ap(0j 0) 2)) (10) 



dt 

dP(0, 1) 
dt 

dP(l, 1) 
dt 



= A 



2A^ 



q 



[P(0,0,2)-P(0,l,2)], (11) 
-[P(l,0,2)-P(l,l,2)]. (12) 



Here q is number of NNs of a vertex, which we as- 
sume to be constant for all vertices in this homoge- 
neous approach. Furthermore, since homogeneity implies 
P(a,a') = P(<r',a), the evolution of all pair is deter- 
mined given Eqs. (|To|)-( 12 ). 



Up to this point, all pair motion equations are exact. 
Now, we perform the pair approximation [38j 



P(er, a', a") 



P{a,a')P{a',a") 
P(a') 



(13) 



which, after applying normalization and symmetry con- 
ditions, leads to the following non-linear equations for 
the stationary state: 



<p _ q p(l - 
A X(q-l)l-p- 

q 



P = i- 



A(g-l) A 



A + l- 



q-i 



(14) 
(15) 



4 



Solving Eqs.(14) and (15) we obtain 



, A - A c (g) 
9 A-A c (g) + 2' 

where the critical point for general q is given by 

2q 



A c (<7) = 



9-1 



(16) 



(17) 



In Sec. [V] we will show that this result fits well the critical 
point obtained from computer simulations of TCP on 
networks. 



B. Heterogeneous mean field theory 

Application of HMF to the TCP in heterogeneous net- 
works follows the steps performed for the CP. In the het- 
erogeneous case, we differentiate the functions p and (f> 
according to degree, defining 4>k and pk as the probability 
that a vertex of degree k is doubly or simply occupied, 
respectively. The rate equations ^ and ^ are thus sim- 
ply modified within the approximation, taking now the 
form 



d<pk 
dt 

dpk 
dt 



= -<j>k + Xkpk&k 

= Afc(l - 2p k - tf>k)6k, 



where 



P(k'\k)J> k , 
k' 



(18) 
(19) 

(20) 



is the probability that a vertex of degree k receives a 
particle from an active NN. As in the CP, the factor l/k' 
accounts for the random choice of the target neighbor 
for offspring creation from a vertex of degree k' . In the 
case of degree uncorrelated networks, we are led to the 
simplified coupling factor 8fc = = ~}2k4>kP{k) / {k) 
A . independent of k. 
In order to solve Eqs. ( 18 ) and fl9| ), we notice that <f>k 
gives the order parameter <fi = J^k P(k)4>k in the APT. 
The density p = ^2 k P{k)pk is a non-critical field coupled 
to the order parameter that goes to a finite value at the 
steady state even at the critical point or in the subcrit- 
ical phase. We thus assume that it will reach its steady 
state in a way faster than <fi. We have checked that this 
assumption is compatible with our numerical simulations 
(see Fig. [2]) . We thus proceed to perform a quasi-static 
approximation 40] , imposing pk — in Eq. ( 19 1 to obtain 



1 



Pk 



(21) 



which, upon substitution into Eq. ( 18 ), leads to the equa- 
tion 



d ^j- ~ -4> k + -fc(l - 4>k)®k, 



(22) 



that holds for long times. 

Equation ( 22 ) is exactly the same HMF rate equation 



obtained for the order parameter of the CP on heteroge- 
neous networks (see Eq. ([!])), with an effective creation 
rate A' = A/2. Therefore, we can exploit all results pre- 
sented in Sec. [TTj For an infinite network size, the TCP 
will undergo a continuous APT at a critical A c = 2, inde- 
pendently of the connectivity pattern of the network. Un- 
der the additional simplification of degree uncorrelated 
networks, we recover the critical behavior <f> ~ (A — X c )^ 
with j3 — max{l/(7 — 2), 1}, while the characteristic time 
diverges as t ~ |A — A c | _1 . At the critical point, the or- 
der parameter decays as (f> ~ t~ s , with S — (3. In finite 
networks, on the other hand, the same mapping to a one- 
step process performed for the CP ensues from Eq. ( [22] ) 
(with the corresponding mapping A' = A/2), and there- 
fore we expect for the TCP a FSS form given by Eqs. p]). 
We finally note that Eqs. (21) and (22) lead to a sta- 



tionary density of singly occupied nodes given by 



Pk 



Xk4 

W) 



o[(ktf>y 



(23) 



Therefore, at criticality about 50% of the nodes be- 
long to the backbone of singly occupied vertices through 
which dynamic spreads. This result is verified with good 



accuracy in our simulation results, presented in Sec. VI 



V. SIMULATION METHODS 
A. The quasi-stationary state method 

In infinite systems, APTs are sharply defined. In finite 
systems, however, the absorbing state can be visited even 
in the supercritical phase due to stochastic fluctuations. 
The numerical exploration of APTs in this case is a subtle 
problem, that must be handled with suitable simulation 
strategies. Here we investigate the stationary properties 
of the TCP using the quasi-stationary (QS) method |32j . 
which has been successfully applied to determine with 
high accuracy the critical properties of the CP in SF 
networks HO]. In the QS method, every time the 
system tries to visit an absorbing configuration, the cur- 
rent configuration is replaced by one active state picked 
up randomly from the system history. The method is 
implemented by storing and constantly updating a list 
containing M active configurations previously visited by 
the dynamics. An update consists in randomly selecting 
a configuration in the list and replacing it by the present 
active configuration with a probability p r At. After a re- 
laxation time t r , the QS probability P n , corresponding to 
states with n occupied vertices, is computed during an 
averaging interval t a . All relevant QS quantities can be 
then determined from P n as, for example, the density of 
active particles <fi — N^ 1 J2 n n Pn an d the characteristic 
lifespan r = 1/P X [19]. 

The networks used for simulations were generated us- 
ing the uncorrelated configuration model (UCM) [41] . 
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with a minimum degree m — 6 and a hard cutoff k c — 
N 1 / 2 , which leads to the absence of the degree-degree 
correlations in large quenched SF networks |36j without 
self nor multiple connections. System sizes with up 10 7 
nodes were used. Simulations were performed with fixed 
parameters t r = 4 x 10 6 , t a = 4 x 10 7 and M = 200. The 
list with the system history must be updated according to 
the frequency that the system tries to visit the absorbing 
state. The larger the lifespan r the slower the list must 
be updated in the stationary state. In our simulations, 
we used p r = 10~ 3 — 10~ 2 , being the larger values used 
for the smaller network sizes. 



B. Critical point determination: The moment ratio 
method 

A key point in the numerical computation of critical 
exponents in any phase transition is the accurate de- 
termination of the critical point. The usual strategy 
to determine critical points in lattice models considers 
the search for pure power laws in plots ef> or t versus 
N [TJ]. This strategy is however not well suited in the 
presence of strong corrections to scaling, as expected for 
the TCP, given its relation with the CP. In Ref. [20], it 
was proposed an alternative method, using the idea of 
size-independent critical moment ratios |42j in the form 



M™ = 



qs mm 1 y 

From the scaling form of P n , Eq. pi, we have 



n , I n 



k ) = Y—f t 

n—l v v v 



(24) 



(25) 



Defining x = n/Vfl with increments Ax = l/VO we find 



N/n 

£ 



N- 



70*0 Az = 



(26) 



where au = L°° x k f(x)dx. Therefore, the asymptotic 
moments ratios at mean field level are size-independent 
and given by M™ s = a n /a q a s . One can thus estimate the 
critical point taking advantage of this fact: Plotting the 
quantities M% s as a function of A for different values of 
N, the different curves will intersect at the critical point 
A = A c , for all different network sizes [42] . Obviously, the 
validity of this scheme is conditioned to the scaling form 
of P n , that holds only for large systems, even at a mean 
field level M- 



of annealed networks (TSJ, H3], in which all edges are 
rewired between any two dynamical steps, while preserv- 
ing the degree and degree correlations of each individual 
node. This procedure destroys all dynamical correlations 
and thus renders exact the predictions of HMF theory 
O [18]. In practice, simulations on annealed networks 
are carried out in uncorrelated networks by selecting at 
random, every time that a node has to choose a nearest 
neighbor, a vertex of degree k in the whole network with 
a probability kP(k)/(k) [18]. In the following, we con- 
sider SF networks with a degree distribution P(k) ~ fc~ 7 
and 7 > 2. 

The HMF theory discussed in section IV B| predicts a 
FSS explicitly depending on the factor g = (k 2 }/(k) 2 
that, for SF networks with 2 < 7 < 3, behaves as g ~ 
const x [l - £ 3 -t + 2f~ 2 ■ • •] A; 3 ^, w here £ = k /k c . 
So, even in annealed networks the FSS is ruled by strong 
corrections, specially for 7 « 3 and 7 w 2 when they 
become logarithmic. However, since HMF theory is exact 
in this case, the critical point must be A c = 2. We have 
checked that the moment ratio method precisely recovers 
this theoretical expectation, with all moment ratios, for 
different values of q and s and different network sizes, 
intersecting at A c = 2 for any value of 7. 

Figure [l] shows the QS density and characteristic lifes- 
pan for critical TCP on annealed networks with degree 
exponents 7 = 2.25 and 7 = 3.25. The dependence of 
these quantities with the system size neatly deviates from 
the pure power law at criticality given by Eq. jHJ). This 
is however, a simple finite size effect; when the full factor 
g is explicitly included, we obtain a perfect agreement 
with HMF exponents </> ~ (gN)~ s " and r ~ (N/g) s <* 
with S v — S a = 1/2. The exponents obtained by means 
of linear regressions are S v = 0.498(3) and 0.502(3) while 
S a = 0.499(2) and 0.498(2) for 7 = 2.25 and 3.25, respec- 
tively. 



VI. TCP ON QUENCHED NETWORKS 

The excellent agreement between HMF and QS simu- 
lations on annealed networks is expected since the cen- 
tral MF hypotheses are fulfilled: absence of dynamical 
correlations and a distance between nodes £ = 10 im- 
plying the exactness of the one-site approximation. In 
real (quenched) networks, whose edges are frozen and 
do not change in time, we usually still have the small- 
world property, but dynamical correlations are unavoid- 
ably present. An additional feature of TCP is that ac- 
tivity spreads through the backbone of occupied vertices. 
So, a necessary condition for the maintenance of activity 
is the existence of a giant component (GC) of connected 
vertices in the backbone [5 . This component is present 
in our simulations as discussed in the end of this Section. 



Check of numerical methods 



In order to check the efficiency of our simulation and 
numerical methods, we have considered the TCP on top 



1 In an annealed substrate a given node interacts with any node 
of the network including itself due to the connectivity rewiring. 
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FIG. 1. Finite size scaling of the critical density (top) and 
characteristic lifespan (bottom) for the TCP on annealed SF 
networks. Symbols represent simulations; dashed lines are the 
scaling as a function of N expected from Eq. (5| for 7 = 2.25; 
full lines are the HMF theory predictions for 7 > 3, Eq. jSJ, 
or for all 7 when taking into account the full g dependence, 
Eq. Q. Filled symbols were shifted upwards to improve vis- 
ibility. 



An important issue in the simulations is that, differ- 
ently from the standard CP on quenched networks |20j . 
the relaxation of TCP in SF networks depends on ini- 
tial conditions. Figure [2] shows the critical relaxation 
in a quenched SF network (7 = 2.75) for two initial 
conditions. A fully active initial condition leads to a 
metastable phase where the system gets stuck into a 
pseudo-stationary state of very low density. After a tran- 
sient, which can be very long for large sizes, the system 
evolves towards the real QS state. However, if we start 
with an initial condition having a small fraction of ac- 
tive vertices and the remaining ones occupied by a single 
particle, the metastable phase disappears and the system 
goes directly to the real QS state. This dependence on 
the initial condition is easily understood by considering 
that, in the early steps of a simulation with a fully active 
initial condition, annihilation events are predominant due 
to the lack of free space to create new particles. These 
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FIG. 2. Critical relaxation of the density of active vertices 
for QS simulation of TCP on quenched SF networks with 
7 = 2.75 and size N = 1.28 x 10 6 Two initial conditions are 
reported: Fully active network (crosses) and 2% of randomly 
distributed active vertices (circles). The dashed line is the 
power law decay <f> ~ t -4 / 3 predicted by HMF theory. Inset: 
Critical relaxation of the density of singly occupied vertices 
in a semi-log scale. The horizontal line indicates the density 
p = 1/2 predicted by the HMF theory, Eq. (23 1. 



annihilation events fragment the backbone into several 
small components and activity cannot spread until the 
backbone is regenerated, after a very large transient. 

We thus choose to start with a 2% initial density of 
active vertices in all simulations reported in the follow- 
ing. Notice that the density of singly occupied nodes 
goes quickly to values very close to 1/2 as shown in the 
inset of Fig. [2] but the final quasi-stationary state takes 
much longer times in consonance with the quasi-static 
approximation used in HMF calculations of Sec. |IVB| 

Critical points were determined using the moments ra- 
tio method described in section [VC| Figure [3] shows third 
and fourth order moments ratios against creation rates 
for networks of different sizes and degree exponents. All 
curves intersect closely at the same point yielding the 
critical points and moments ratio values shown in Ta- 
ble [T] These moments ratio values are noticeably in good 
agreement with those obtained for CP simulations on the 
same kind of networks [20]. Notice also that the critical 
points are larger than the HMF prediction A c = 2, in- 
dependently of 7. A satisfactory agreement between the 
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Mi 2 


2.50 
2.75 
3.25 


2.1280(5) 
2.1661(2) 
2.2229(1) 


1.87(2) 
1.790(6) 
1.712(8) 


2.65(2) 
2.43(3) 
2.27(1) 


4.73(5) 
4.15(9) 
3.71(4) 



TABLE I. Critical points and moment ratios for TCP on 
quenched SF networks. 
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FIG. 3. Moment ratios against creation rate of the TCP on quenched SF networks with different values of 7, for several system 
sizes. For the sake of comparison, all plots have the same horizontal range. 
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FIG. 4. Comparison between critical points obtained in a ho- 
mogeneous pair approximation (symbols) and QS simulations 
(lines) of the TCP on quenched SF networks. HPA results 
correspond to Eq. ( |17[ ) with the coordination q replaced by 
the average degree of the network. 



homogeneous pair approximation value for the critical 
point, Eq. ( |17[ ), and the QS simulations can be however 
observed, see Fig. |4j again in conformity with the results 
for the CP on quenched networks [301 13H] • 

Now in possession of accurate estimates of the crit- 
ical points, we can check the validity of HMF critical 
exponents. Figure [5] shows the QS density and charac- 
teristic time for 7 = 2.50 around the estimated critical 
point. Performing power law regression of these quan- 
tities as a function of N, a very good agreement with 
HMF exponents in Eq. ^ is found, as shown in Ta- 



ble [TTJ Using the full anomalous FSS predicted by HMF 
theory, performing power law regressions of the form 
~ (Ng)~ Sv and r ~ (N/g) Sa , all numerically obtained 
exponents are in very good agreement with the HMF re- 
sult S a = S v = 1/2. The agreement of simulations with 
HMF theory extends to all the values of 7 considered, as 
shown in Table ITT1 

We have finally investigated the backbone of the net- 
work used for the spreading of activity, which is defined as 
the largest connected cluster (giant component) of ver- 
tices having at least one particle. Figure [6] shows the 
average fraction of vertices belonging to the backbone in 
the critical TCP against network size. The size of the GC 
corresponding to a random dilution of 50% is also shown 
for the sake of comparison. It is well known that SF 
networks are resilient to random removal of nodes, such 
that a giant component persists even if a finite fraction 
of the nodes is removed [TTJ [12] . A random dilution of 
50% of the nodes in UCM networks leads to GCs contain- 
ing almost the totality of the non-diluted nodes (Fig. 6|. 
The critical TCP dynamics generates a backbone with a 
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VHMF 


s v 


a 


&HMF 


Sol 


2.50 
2.75 
3.25 


0.62(1) 
0.57(1) 
0.50(1) 


0.625 
0.5625 
0.500 


0.51(1) 
0.496(8) 
0.48(1) 


0.37(1) 
0.42(1) 
0.49(2) 


0.375 
0.4375 
0.500 


0.48(1) 
0.50(1) 
0.51(2) 



TABLE II. Critical exponents for TCP on quenched SF net- 
works. HMF exponents Ohmf and &hmf, as given by Eq. jHJ, 
are included for comparison. 
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FIG. 5. Finite size scaling of the QS density (j> and characteris- 
tic time r around the critical point for quenched SF networks 
with 7 = 2.50. Solid lines are linear regressions used to es- 
timate the scaling exponents. Open symbols refer to plots 
against size TV while filled symbols to the plots against size 
re-scaled by the factor g = (k 2 )/(k) 2 . 
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FIG. 6. Average size of the giant component against network 
size obtained with the TCP dynamics or a random dilution of 
50% of the vertices. The random dilution for other 7 values 
have the same behavior of 7 = 2.50. 



similar fraction of vertices. The second largest connected 
cluster contains just a few nodes (typically 3 or 4 for a 
network size N ~ 10 7 ) and, therefore, the activity in 
regions disconnected from the backbone is negligible. 

Noticeably, while random dilution and critical TCP 
on homogeneous and large networks (7 = 3.25) generate 
GCs that do not vary with network size, in the hetero- 
geneous case we have a slow (sub-logarithmic) increase, 
which is more apparent the smaller the value of 7. This 
dependence of the backbone size with IV is however so 
slow (e.g., a relative variation of around 0.5% when size 
increases from 10 3 to 10 7 for 7 = 2.50) that its effect is 



essentially negligible, as can be seen from the excellent 
agreement with HMF theory shown by numerical simu- 
lations. 



VII. DISCUSSION 

In the present work, we have studied the absorbing 
state phase transition of a threshold contact process 
(TCP) in networks. The model, possessing infinitely 
many absorbing configurations in the thermodynamic 
limit, is analogous to the pair contact process [24] with a 
relaxation of the fermionic constraint: up to two particles 
can occupy a vertex but only the particles in doubly oc- 
cupied vertices react (annihilate or create a new particle 
in a nearest neighbor). 

Applying the standard heterogeneous mean-field the- 
ory and extensive large-scale numerical simulations based 
on the quasi-stationary state method, wc have shown 
that the critical behavior of the TCP is exactly the same 
as the standard CP, and therefore that both models be- 
long to the so-called CP-HMF universality class in net- 
works, which is defined, at the FSS level, by the set of 
exponents given by Eq. (|5j). At this respect, the apparent 
non mean-field exponents reported in Ref. [2 9) for the re- 
lated pair contact process, a fermionic counterpart of the 
TCP, can probably be attributed to the networks sizes 
considered in that work (N < 10 4 ), much smaller than 
the largest system sizes (N ~ 10 7 ) that we attain here. 

The CP-HMF universality class studied here can be 
conjectured to include other CP-like models with a fi- 
nite or infinite number of absorbing states. However, its 
status is not as robust as in lattices, since in networks 
microscopic details of the models can be enhanced due 
to heterogeneity. For example, the susceptible-infected- 
susceptible (SIS) epidemic model 44 is a variation of the 
CP in which a vertex produces activity in all its neigh- 
bors with the same strength, while in the CP the activity 
production of a vertex is equally divided among all neigh- 
bors (factor 1/k' in Eq. (pi)])). According to the Janssen- 
Grassbergcr conjecture [45, 0Hj , the universality class of 
CP and SIS should be the same, and in fact both mod- 
els belong to the DP universality class in homogeneous 
lattices. In complex networks, however, heterogeneity 
strongly affects both models and it renders the SIS with 
critical exponents at the HMF level which are different 
from the CP [H]. Moreover, the SIS can be shown to ex- 
hibit a vanishing critical point in the thermodynamical 
limit if the maximal connectivity is unbounded |47[ 148] . 

We can thus conjecture a splitting of universality 
classes, comprising at least a CP-HMF and a SIS-HMF 
class, which are independent of the number of absorb- 
ing states (as proved in the SIS case by the study of the 
forest-fire model [28]) but depends on the way in which 
activity spreads over nearest neighbors. Further work 
should be devoted to the study of this universality class 
splitting in complex networks. 
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